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Abstract 

A hierarchy of mathematical models describing viscosity-stratified flow in a Hele- 
Shaw cell is constructed. Numerical modelling of jet flow and development of vis¬ 
cous fingers with the influence of inertia and friction is carried out. One-dimensional 
multi-layer flows are studied. In the framework of three-layer flow the interpretation 
of the Saffman-Taylor instability is given. A kinematic-wave model of viscous fin¬ 
gering taking into account friction between the fluid layers is proposed. Comparison 
with calculations on the basis of two-dimensional equations shows that this model 
allows to determine the velocity of propagation and the thickness of the viscous 
fingers. 
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1 Introduction 

A displacement process involving two fluids is often unstable when the displacing fluid 
has larger viscosity than the displaced one. The resulting instability developing at the 
interface between two fluids is known as viscous hngering [201 [ig. This instability has 
received much attention as an archetype of pattern-formation problems and as a limiting 
factor in the recovery of crude oil. Classical mathematical model describing a Newtonian 
how displacement in a Hele-Shaw cell and development of the Sahman-Taylor instability 
consists of the continuity equation, Darcy’s law and a convection-dihusion equation for 
the concentration of the displacing huid [221 13]. The inertia of huid may be important 
for high huger velocities. This leads to the necessity to use more complex nonlinear 
equations of huid motion HD]. In the framework of these models instability caused by 
diherent velocities of layers movement can be considered. There are several theoretical 
and experimental studies on the role of inertia in immiscible laiH] and miscible 1271 
displacements. The results reveal that inertia tends to damping viscous hngering. In 
recent publications mm the ehects of the thickness variation of a Hele-Shaw cell and 
elasticity of the walls on the process of viscous hngering have been studied. Diherent 
types of instability in viscosity-stratihed how have been discussed in HU. 

A number of theoretical, numerical, and experimental works devoted to understanding 
of various aspects of the instabilities is available in the literature [TBl [23112T] . Accurate 
numerical solution is costly at high Peclet numbers and it is difhcult to reproduce the 
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detailed fingering pattern. Nnmerical simulations [281 ES] show that it is simpler to 
describe the concentration of solute averaged across the fingers. In such cases, the mixing 
zone is an important feature to determine the extent of mixing. Despite the considerable 
work done, the spreading and growth of the mixing zone is an important question that still 
remains unresolved. Several empirical models are available for the evaluation of mixing 
zones in unstable, miscible displacements. Two empirical models have been suggested 
by Koval [I3] and Todd and Longstaff 121 to give a basis for computation of miscible 
displacement. Both models suffer from adoption of empiricism in which the principal 
parameters involved have little or indirect physical significance. Further development of 
averaged models of fingers formation is represented in [HlEHlIl]. All these models are based 
on the hypothesis of pressure equalization in the transverse direction to the main flow, 
as well as an empirical information about the displacing and displaced fluids distribution 
in the region of intensive viscous fingering. Let us recall that 2D numerical solution with 
high resolution is hard to construct, that is why ID models play an important role in 
some cases. For instance, these models are very useful in calculating of fracturing when 
it is necessary to solve the equations of the crack opening and the fluid motion in the 
fracture simultaneously [I]. 

The aim of the present paper is to derive a hierarchy of mathematical models describing 
viscosity-stratified flow and spreading and growth of the mixing zone in a Hele-Shaw cell. 
In Section 2 we propose 2D nonlinear hyperbolic system of balance laws. In contrast with 
widespread model of flow displacement in a Hele-Shaw cell we apply nonlinear momentum 
equations and take into account compressibility of the fluid. At the same time diffusion 
coefficient is neglected that corresponds to the large Peclet number limit. As we show in 
Section 3 by numerical calculations, this model is suitable for describing of jet flow and 
propagation of viscous fingers in a Hele-Shaw cell. We also point out that for the process 
of unidirectional displacement the pressure variation in the transverse direction is small. 
This observation makes it possible to use long-wave approximation and construct a class 
of layered flows described by a system of one-dimensional evolution equations. Based on 
the various simplifications of the momentum equation (linearisation, lubrication theory) 
a hierarchy of ID mathematical models is constructed in Section 4. The equations of 
a three-layer flow are studied and numerical computations of the formation of viscous 
fingers are performed in Section 5. We show that three-layer stationary flow is correctly 
described in the framework of simplified model. Nevertheless the growth rate of viscous 
fingers is significantly higher than it is observed experimentally. In section 6 we propose 
ID kinematic-wave model of viscosity-stratified flow taking into account friction between 
the fluid layers. The velocity of propagation and the thickness of the viscous finger in the 
framework of this kinematic-wave model coincide with the corresponding calculations on 
the basis of the 2D equations. It gives possibility to predict the parameters of viscous 
fingers without time-consuming calculations. We also show that the proposed ID model 
is in good agreement with the well-known Koval model. 
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Figure 1: Hele-Shaw cell geometry. 


2 Mathematical model 

A Newtonian weakly-compressible flow displacement in a Hele-Shaw cell (the area between 
two parallel plates separated by a small gap of constant thickness b in the z direction, see 
Fig. [T]) is described by the equations 

(pv)i + div(pv ® V - P) = 0, pt + div(pv) = 0, = 0- (1) 

Here v = (n, v, w) is the fluid velocity, p is the density, and P is the stress tensor. This 
tensor can be taken in the form P = — (p -|- |/idivv)l -|- (Vv -|- (Vv)*)/i, where p is the 
pressure, and p is the viscosity. To describe the process of displacement involving two 
fluids of different constant viscosities we should take into account that viscosity p depends 
on the concentration of solvent c. This function is scaled such that it is equal to unity 
in the displaced fluid (p = P 2 ) and zero in the displacing one (p = pi). Following [3] 
we assume a monotonic relationship between the viscosity and the concentration in the 
form p(c) = p}“'^p 2 . The concentration c satishes to the transport equation (diffusion is 
neglected that corresponds to the large Peclet number limit) 

Q-|-v-Vc = 0. (2) 

We assume now that the velocity held can be represented as 

M = ^(1 - (y) )M'(t,x,p), V = ^(l - w = 0. (3) 

It provides the fulhlment of no-slip conditions on the cell walls 2 ; = ±6/2. We also suppose 
that the functions p, p, and c do not depend on 2 ;. Let us note that in the calculation of 
divP the following terms and {p^y)y can be omitted as they are negligible 

compared to the derivatives with respect to z. Further, averaging Eqs. ([I]) and ([2]) through 
the gap we obtain (primes are omitted) 

{pu)t + {Ppu^ + p)x + {/3puv)y = -pu, 

{pv)t + {I3puv)x ± (/3pu^ ± p)y = -pv, (4) 

Pt + {up)x + {vp)y = 0, (cp)i ± {ucp)x ± {vcp)y = 0. 
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Here and below /i denotes the viscosity of the fluid divided by the permeability 6^/12 
(further in the text we will call it simply “viscosity”); coefficient /3 is equal 6/5 (this 
factor comes from integration of v 0 v in the form ([3]) with respect to 2 ; [inii). 

In order to close model dl]) we should specify either the equation of state p = p{p) 
(barotropic fluid) or dependence p = p{c) (incompressible fluid). A weak compressibility 
of the fluid given by the equation of state p = p{p) provides hyperbolicity of the model. 
This dependence can be considered as a regularization of equations describing the flow 
of an incompressible fluid in a Hele-Shaw cell. On the other hand in applications the 
property of hyperbolicity of equations can be related to the presence of gas cavities in the 
porous medium as well as with the elasticity of the channel walls (for instance in PKN 
model [1] the pressure p depends on the channel thickness b). In any case if the condition 
p'{p) holds then results weakly depend on the choice of p = p{p). Therefore for 
the numerical simulation of 2D flows we assume 

p{p) = aVV2 (a^ = cl/po) (5) 

where the constants po and cq specify characteristic density and speed of sound in the 
fluid. 

To hnd the characteristics of system dH), ([5]) we write it in the vector form 

Ui + AU^ + BUy = F 

where U = {u,v,p,c)'^ is the vector of dependent variables; F = {—pu/p,—pv/p,0,0)'^ 
is the right-hand side; A and B are 4x4 matrices. Let ^ = (G;‘^ 2 , ■Cs) be the normal 
vector to the characteristics; I is the identity matrix. Then the characteristic matrix 
C = ^il -I- ,^2 A -I- ^sB of system (jl]), ([5]) has the form 

^xi + (/^-i)w6 (/5-i)w6 ((^-i)(x2-6)w + p'(p)6)p“^ 0^ 

^ (/ 5 - 1 X 2 xi + (/^-iX3 ((/^ - i)(x 2 -6)^^+ p'(p)6)p”^ 0 

X2 0 

V 0 0 0 X2j 

Here Xi = ■Ci + X 2 = ii + u ^2 + A simple but cumbersome calculation 

yields the following expression for det C(^): 

det C(|) = ((^1 + 2 / 3 (m ^2 + < 3)6 + f3{u^2 + < 3 )^) - (^^ + v'^)p'{p))xiX 2 - 

We specify the characteristic surface by the equation W{t,x,y) = 0. Then to obtain 
the differential equations of the characteristics we should replace the vector (■Ci,'C 251 ^ 3 ) 
previous equation by the vector (Wt, HA, Wy) and equate to zero det C. As a result we 
obtain two families of contact characteristics 

Wt + uW^ + vWy = 0, Wt + + PvWy = 0 

and two additional characteristic families 

Wt + PuW^ + l3vWy = - l){uW, + vWyY + (WJ + w^)p'{p). 
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If the inequalities /3 > 1 and p'{p) > 0 hold this system of equations is hyperbolic. Note 
that in the case of /3 = 1, /i = 0, and c = const system (jl]), ([5]) coincides with the 
well-known shallow water equations. 

3 Modelling of viscous fingering and jet flows 

Below we present the results of numerical calculations of the viscous hngers and jet streams 
on the basis of hyperbolic model (jl]), (|1]). Originally mathematical description of the 
Saffman-Taylor instability was given in the framework of the Darcy’s law and the mass 
conservation equation [201122]. The inertia of the fluid may be signihcant for high huger 
velocities. In [7] simulation of viscous hngers was performed using non-linear equations 
of an incompressible huid. We show that Eqs. (|1|), ([5]) taking into account the forces of 
inertia and compressibility of the huid could be also used for description of this instability. 

To solve diherential balance laws (|4]), (|5]) numerically one can apply methods based 
on various modihcations of Godunov’s scheme. In this work we implement the robust 
and stable Nessyahu-Tadmor second-order central scheme ng. In every test we assume 
that on the boundaries ?/ = 0 and y = H the impermeability condition n = 0 is fulhlled. 
The size of the computational domain is L = 100, H = 50; the resolution of the problem 
on the X and y axes are 300 and 150 nodes correspondingly (uniform grid). We assume 
that po = 1, Co = 150 and (3 = 6/5 (for the third test (3 = 1). The values of the 
variables are considered as dimensionless. Below we present calculations showing the 
possibility of modelling the evolution of perturbations caused by the Kelvin-Helmholtz 
and/or Sahman-Taylor instabilities on the basis of the hyperbolic model l]T|). fl51). 

3.1 Test 1. Jet flow 

Let at the initial time t = 0 the Hele-Shaw cell be occupied by a quiescent fluid having 
density p = I and viscosity pi = 0.1. Through the left central cross-section of width 
H/10 fluid of viscosity p 2 = 0.4 is injected with velocity U 2 = 24; through the rest part 
of the left boundary fluid of viscosity pi = 0.1 is entered with velocity Ui = 4. On 
the right boundary of the domain the condition of constant pressure is valid. For more 
intensive development of the perturbations at each time step we slightly disturb boundary 
conditions at a; = 0. Namely, the cross section, through which fluid “2” is injected, is 
randomly shifted up/down from its initial position on hxed distance Ay which is equal 
to the grid spacing with some positive integer factor k. We call it “random shake” and 
take here k = 1. The function c for values of the concentration is presented in Fig. [2] 
at t = 10. Vortices are formed at the interface of the layers due to Kelvin-Helmholtz 
instability. Increase of the both values pi and p 2 suppresses this instability. 

3.2 Test 2. Viscosity-stratifled flow 

Viscosity-stratihed multi-layer flow is shown in Fig. (Hjat t = 15. Initially the flow region 
is hlled by a quiescent fluid (pi = 1, p = 1). Liquid is pumped through the left boundary 
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Time 10 


Time 15 





Figure 2: Kelvin-Helmholtz instability in 
jet flow of viscosity-stratified fluid injected 
at X = 0 (/ii = 0.1, Ui = 4, and /i 2 = 0.4, 
U 2 = 24). 



Figure 3: Modelling of multi-layer jet flow. 
Fluid of viscosity /ii = 1 and fi 2 = ^ (for 
odd and even layers) is injected at x = 0 
with velocities Ui = 21, and U 2 = 7. 


divided into layers of height H/12] velocity and viscosity for odd and even layers are 
f /2 = 7, /i 2 = 3 and Ui = 21, /xi = 1, correspondingly. As in the previous example 
“random shake” of the jets is used. The results of the calculations show that multi¬ 
layer flow without mixing is realized for a wide range of parameters (Kelvin-Helmholtz 
instability at the interfaces between the layers occurs if viscosity decreases more than in 
hve times). We note that the condition on the left boundary Ufi = const corresponds to 
a class of exact solutions of equations ([1]) for an incompressible fluid: u = U{y), v = 0, 
p = const, p = —ax, p = a/U{y). The stability analysis of this class of flows was carried 
out in [B]. 

3.3 Test 3. Viscous fingering 

The following example illustrates the formation of viscous hngers. Let the displacing 
phase injected at a constant velocity U be referred to with index 1 and the displaced one 
with index 2. At t = 0 fluid “1” is located in the domain x < Xq = more viscous 
fluid “2” — in the domain x > Xq. For convenience we use the coordinate system moving 
with velocity U and assume /S = 1 (with /3 = 6/5 results are similar). Let us perturb 
the initial interface x = xq as follows: x = xq -l- A?/ cosibuy/H) (here A?/ = 2/3). We 
note that the physical effect of the instability can be obtained numerically if the initial 
perturbation is not less than the grid resolution. At t = 0 we choose piece-wise linear 
pressure distribution (p^. = —piU for x < Xq and = —P 2 U for x > Xq, on the right 
boundary p is equal to CqPo/ 2). Initial density of the fluid is determined using formula (jS]). 
At the boundaries the impermeability condition is fulhlled. 

The calculations are performed for the following parameters: U = 3, pi = 1, and 
P 2 = 5. In the evolution process of the flow viscous hngers are formed (Fig. 01 left). The 
number of hngers is determined by the initial perturbation. Displacing huid penetrates 
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Figure 4: Formation of viscous fingers in weakly-compressible Hele-Shaw flow for U = 3, 
Hi = 1, and /i 2 = 5: the concentration c (on the left) and the density p (on the right). 


more rapidly into displaced one (fingers are not symmetrical with respect to the initial 
interface). Fig. 0] (right) shows the distribution of the density at f = 30. Calculation 
with better resolution leads to the same result. As we can see the density changes less 
than 0.5% in comparison with the initial one (to reduce the compressibility we should 
increase the speed of sound cq but this slows down the calculations since the time step is 
determined by the Courant number). Note that the density (pressure) varies slightly with 
respect to y. This allows one to use approximate model, where the second momentum 
equation is replaced by Py = 0. 

4 Layered flows 

We consider an incompressible {p = 1) Hele-Shaw flow in a two-dimensional domain 
of rectangular geometry (dimensions L and H, respectively) governed by Eqs. (jl]). It is 
assumed that the flow is essentially parallel and the pressure gradient in the flow direction 
X being independent of the transverse coordinate y to leading-order in the variable e = 
H/L. This regime is termed as the state of transverse flow equilibrium [26] . Such flows can 
be also considered in the framework of long wave approximation mu- Let us perform 
the following scaling in Eqs. (jl]) 

t —)■ X e~^x^ V ^ ev, /i —)■ ep. 

Then we neglect terms of order <^1. As a result we obtain the approximate model 

Ut + f^uUx + (3vUy + Px = —pu, Py = 0 , 

Ux + Vy = 0, Ct + UCx + VCy = 0, (6) 

v\ „ = 0 , v\ jj- = 0 , 

\y=0 ’ \y=H ’ 

where the pressure does not depend on the variable y. We also suppose that on the 
boundaries y = 0 and y = H the impermeability condition is fulfilled. 
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Let us consider the class of viscosity-stratified flows 


u = Ui{t, x), c = Ci = const, y G (l/i-i, ?/*) 


(0 = I/O < yi{t,x) < ... <yN = H). In this case Eqs. ([6]) take the form 


Uit +(3uiUi^+p^ =-yiUi, hit + {uihi)^ = Q, = 

N N 

^ ^ hi H, ^ ^ Uihi Q. 
i=l i=l 


( 7 ) 


Here hi{t,x) = yi{t,x) — yi-i{t,x) is the depth of i—th liquid layer of viscosity /ij having 
velocity Ui{t, x); and Q is the total flow rate through the cell. Upon derivation of Eqs. ([7]) 
the kinematic condition at the layers interface is used. 

Introducing new unknown variables s* = Ui — allows to transform Eqs. ([7j) to the 
evolution system of 2(iV — 1) equations 


^it “1“ /^(('5i/2 “1“ (hAf 

hit “1“ ”1“ '^N^hijx 0, (i 1, ..., N 1) 

where 

N-l ^ N -1 

hj\i 1 ^ ^ hi., U]s[ ^ ^ Sihj^. 

i=l i=l 

Further we assume that Q = const. 

We also use the following simplihed versions of governing Eqs. ([7]). The hrst one 
consists in the linearisation of the momentum equations: 

T “t“ Px f; •••; -^) 

(here U = Q/H is the average velocity). The second simplihcation is based on the Darcy 
law: 

Px = -yiUi {i = l,...,N). 

The remaining equations of system ([7]) do not vary. 

In some cases it is convenient to use a moving coordinate system x' = x — Ut, u[ = 
Ui — U. Then Eqs. ([7]) take the form (primes are omitted) 


Uit + (/^Mi + (/5 - h)U)uix +Px = -fJ^iiui + U), 

N N 

hit T {^Uihijx 0? ^ ^ hi H, ^ ^ Uihi 0. 

i=\ i=l 

Further we show that in the framework of three-layer and two-layer regimes of flow it is 
possible to give an interpretation of the Saffman-Taylor instability as well as to describe 
the initial stage of viscous hngering. 


5 Three-layer flow 

We introduce the following notation for the layer velocities and depths 
M = Ml, V = U 2 , w = M3; h = hi, r] = h 2 , ( = ^3- 
We also assume that hf = 1, Q = 1, = pa = 1, and p 2 = h 7^ 1- 


5.1 Stationary solutions 

Here we construct a steady-state solution of model ([7]) for three-layer flow. Integration of 
the equations of conservation of mass in system ([7]) allows to express the depths of the 
layers 

h = Qi/u, r] = Q 2 /v, ( = Q^/w. (9) 

Here Qi is the flow rate in the i-th layer (Qi + Q2 + Q3 = !)• Due to the unit depth we 
obtain the velocity in the intermediate layer 


V = (p{u, w) 


Q2 A _. Qi Qs 

A ’ “ M m; ■ 


Eliminating pressure p from the equations 


(3uu' + p' = —M, (3vv' + p' = —pv, (3ww' + p' = —w 


(here the prime denotes the derivative with respect to x) reduces the problem to the 
solution of the autonomous system of ordinary differential equations 

du {u — pip)w — {u — w)(p(pu, dw {u — pip)u + {u — w){(p(pu — u) 
dx {{(fpu - u)w + uipipw)(3 ’ dx {{ipipu - u)w + uip(pw)(3 

A hxed point of system flTUD is determined from the relations u = w = pp: 


= 1 + {p - 1)Q2- 

Linearisation of Eqs. CO on the solution u = u^, w = and computation of the 
eigenvalues of the corresponding matrix show that the hxed point is a stable node. The 
integral curves in the phase plane (m, w) in the neighbourhood of the hxed point are shown 
in Fig. Ofor Qi = 0.4, Q 2 = Qs = 0.3, p = 2, and fd = 6/5. 

Fig. E] shows a comparison of the numerical results obtained on the basis of two- 
dimensional hyperbolic equations (jl]), ([5]) and multilayer model ([7]) reduced to dynamical 
system flTOll in the case of three-layer stationary how. Solid white lines indicate the 
layers depths y = h{x) and y = h{x) + p{x) obtained by solving equation ffTOj) and using 
relations (E]); dotted lines correspond to the hxed point (h* = Qi/m*, = 1 —Qs/tc*). 

Here we take the following values of layers depths = pq = 0.2, (q = 0.6 at x = 0. As 
before we choose Qi = 0.4, Q 2 = Qs = 0.3, p = 2, and {3 = 6/5. The hgure shows that 
the solution reaches an equilibrium state for x > 5. 
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Figure 5: The integral curves of ODE ffTOj) 
in the phase plane (a stable node) obtained 
for Qi = 0.4, Q 2 = Qs = 0.3, and fi = 2. 



1 2 3 4 5 6 7 


JC 

Figure 6: Comparison of the results for 2D 
Eqs. (jl]), ([S]) (concentration c is shown in 
two colours) and Eqs. fITUD for three-layer 
flow (solid white lines). 


To carry out the calculation on the basis of 2D equations (j4]), ([5]) the following initial 
data are used. The flow domain in the ^/-direction is divided into three layers of width 
ho, rjo, and Co- The fluid of density p = 1 at f = 0 moves in these layers in the x-direction 
with constant velocities u = Qi/ho, v = Q2/vo^ and w = Q3/C0 respectively. We also 
suppose that p = 1 in the layers of width ho and Co; in the middle layer of width rjo we 
choose p = 2. The same data are taken as the boundary conditions at x = 0; on the right 
boundary (x = 8) the condition of constant pressure is prescribed; on the walls y = 0 and 
y = I the condition of impermeability is fulhlled. The resolution of the problem on the x 
and y axes are 300 and 60 nodes correspondingly (uniform grid). In order to visualize the 
flow of fluids with different viscosities the concentration c is used. This value is presented 
in Fig. |6]at t = 25. More viscous fluid in the middle layer is shown in brown (c > 0.35) 
and less viscous one is shown in blue (c < 0.35). 

5.2 Non-stationary solutions 

Let us consider a three-layer flow governing by Eqs. ([7]) wherein the momentum equations 
are replaced by linear Darcy laws Px = —piUi. Taking into account assumptions above 
and notations we have 

u = w = p/d, v = l/d, d = {1 — p)ri + p. (11) 

In this case the depths of the layers h and rj are found from the system of equations 

ht + {u{r])h)x = 0, pt + {v{r])p)x = 0. (12) 

It is easy to check that this system is hyperbolic and its characteristic velocities are 

\l=r]v'{7]) +V{7]), \2 = u{7]). 
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The first family of characteristics is genuinely nonlinear whereas the second one is linearly 
degenerate [I9]. In terms of the Riemann invariants rj and r = h/{l — rj) Eqs. (1121) take 
the form 

Vt +>^i{v)Vx = 0, rt + X2{r])rx = 0- 

In the case 0 < /r < 1 we construct a centred simple wave solution dehned by the 
relations 

r = ho = const, Ai(? 7 ) = ^, ^ = {x — xo)/t 

(note that the ansatz 7] = const leads to a constant solution). The layer depths are 


Formulae (IT^ give the solution of Eqs. (IT^ with discontinuous initial data 




t=o 


(0, 1), X <xo 

{ho, 0), X > xo- 


(13) 


(14) 


Prohles of a viscous huger y = h and y = h + rj given by flT3|) are shown in Fig. [7] for 
ho = 0.6 and various values of /r < 1. 

Let us construct a solution of Cauchy problem flTT]) . (ITT)) for /i > 1. At initial time 
the velocities and the layers depths are u~ = y, u'^ = 1, = 0, = ho] v~ = 1, 

= l/y, r]~ = 1, ? 7 + = 0. It is easy to verify that these values satisfy the Hugoniot 
conditions 

[{u — D)h] = 0, [(n — D)r]] = 0 

derived from Eqs. (IT^ as well as the stability conditions [19] if the shock front moves 
with average how velocity D = U = 1. 

It is interested to note that in the class of simple wave solutions (r = const) system (1121) 
reduces to the naive Koval model HISS]. In fact, taking into account representation (fTTD 
the second equation in (fT2l) can be rewritten as 

dc d f Me \ /icrN 


where c stands for y and M = 1/y. It is known that the growth rate of viscous hngers 
in the framework of the naive Koval model is signihcantly higher than it is observed 
experimentally. 

To derive another model of a three-layer how in a moving coordinate system we use 
linearisation of momentum equations in (151) 


uu + jUui,x+Px = -hi{ui + U), (7 = /3-1). 


Eliminating the pressure p leads to the system of evolution equations 


Sit + ISlx — —Si, 

S2t + 7S2x = -hS2 + {l-p){l-w), (16) 

ht + (uh)^ = 0 , r]t + {vr])x = 0 
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Figure 7: Layers thickness y = h and y = 
h + r] obtained by self-similar solution flT3|l : 
1 — yU = 0.3, 2 — /i = 0.5, 3 — yU = 0.7. 



Figure 8: Profiles ol y = h and y = h + r] 
(solution of flTB]) is given by solid curves, 
dot-dash corresponds to ffT3|) : dashed line 
presents the initial data): 1 — t = 2, 2 — 
t = 5. 


where 


U = Si + W, V = S2 + W, W = —hSi — 7jS2- 


Let us rewrite Eqs. (IT 6 |) in the form Ut -|- AUj, = F. Here U = (si, S 2 , h, is the 
unknown vector, A is the Jacobi matrix, and F is the right-hand part. The eigenvalues 
of matrix A are 

Ai = ysi, A2 = 7S2, 

A34 = 2“^((1 - 3 h)si + (1 - 3 ? 7 )s 2 ± \/m) 

where m = (s 2 — (1 — h)si)^ -|- -|- 2((1 -|- h)si — S 2 )r}S 2 - Conditions 0 < h < 1 and 

0 < 7 < 1 provide hyperbolicity of system flTOj) since m > 0. Introducing the parabolic 
with respect to r = Si/S 2 function / = mjs\ it is easy to check validation of inequalities 
/"(r) > 0 and /(r*) > 0 (here r* is a minimum point of /(t)). It means that m > 0 and, 
consequently, characteristic velocities A* are real. 

Further we construct the numerical solution of Eqs. ffTHj) with initial data 


[ (0,1 - /i,0, 1), X < xo 

(si, S 2 , h,7) , n = < 

I (0,(1 -yu)/y(i, ho, 0), X > Xq. 

Note that this formulation corresponds to Canchy problem flT^ . ffT4)) . Calculations on 
the basis of model ffTHj) are carried out using Nessyahu-Tadmor scheme na. The results 
of computations are shown in Fig. [8] at various moments of time (solid curves). As we 
can see with increasing of time the solution tends to the self-similar regime (dot-dashed 
curves obtained by using formulae flT3l) L Moreover tip of the viscous finger propagates 
with the same velocity in the framework of models ffTHj) and flT^ . 
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6 Kinematic-wave model 


In the previous section it is shown that three-layer flow governed by the simplifled ID 
model fll2p (or correctly describes the well-known fact that the fluid interface is 

unstable if less viscous fluid displaces more viscous one and stable otherwise. However, 
the velocity of the viscous Angers propagation for these equations is the same as for 
the naive Koval model which vastly over-predicts the rate of the Angers grow [U [26]. 
A number of empirical models is proposed to reconcile this behaviour [T21 El El]- For 
example, Koval [13] postulates empirically that flT3]) is valid if the viscosity ratio M is 
replaced by an effective viscosity ratio Me, where 

Me = (M'/^Ce + 1 - Ce)^ Ce = 0.22. (17) 

Despite the fact that this model appears to work in practice [I5] a theoretical justification 
is yet to be obtained. 

As it can be seen from Fig. [7] and [SI in the framework of model 012!) or fllOp . the 
tip of the viscous Anger becomes infinitely thin with time. However, in experiments the 
formation of Angers of finite thickness is observed. The friction between the fluid layers 
is one of the possible mechanisms that prevent thinning of the viscous Anger tip. Below 
we obtain modification of the above considered ID models, which takes into account 
friction between the fluid layers. This model is new and correctly describes growth of 
the viscous Angers. It is proved by comparing with the Koval model (Eq. flTbll where an 
effective viscosity ratio Mg is used) as well as by numerical calculations on the basis of 
2D Eqs. dH), ([3]). 

6.1 Friction between the layers 

Let us consider a two-layer flow governed by the following system of equations 

Ut + {3uux + Px = —piu — Kfl{u — v)h~^, 

Vt -F f3vVx +Px = -P 2 V K,fl{u - v)p~^, 

(18) 

ht + {uh)x = 0, Pt + {pv)x = 0, 
h + p = 1, uh + vp = 1. 

Here we use notations of the previous section for the layers velocities and depths; constant 
K > 0 is an empirical parameter, and p = y/piP 2 - This system differs from Eqs. ([7P (in the 
case iV = 2) by the presence of additional terms with factor pn which expresses friction 
on the layers interface. Note that this effect becomes significant if the thickness of one of 
the layers tends to zero. 

Further we use simplification of the momentum equations in flTSp based on the Darcy 
law: 

Px = —piU — Kjl{u — v)h~^, Px = —p 2 V -|- KfL{u — v)p~^. 

Taking into account that 

p = 1 — h, n = (l — uh)p~^, piu + Kfl{u — v)h~^ = p 2 V — Kfl{u — v)p~^ 
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ffT^ with M = 10, K = 0.45 and its convex 

hull. 



Figure 10: Comparison between the Koval 
model (solid curve) and model flT^ with 
K = 0.45 (dashed) for M = 4. 


we obtain the following kinematic-wave model 


ht + (<h(h))a, = 0, ^ = uh 


(k\/M+ (1 - h)hM)h 

(1 - h){l - (1 - M)h)h ’ 


(19) 


where M = Typical graph of the flux ^{h) is given in Fig. [9] (curve 1) for M = 10, 

and K = 0.45. Obviously *h(h) is a non-convex monotonic function such that *F(0) = 0, 
and <h(l) = 1. 

We construct self-similar solution of Eq. flT^ with initial data h(0, x) = 1 for a; < xq 
and h{0,x) = 0 for x > Xq. On the interval h G (0,1) function ^{h) has two points of 
inflection. Thus it is necessary to construct convex hull of the flux 4)(h). For this we draw 
tangents to 4’(h) from the origin and from the point (1,1) (lines 2 ans 3 in Fig. IH]). Let 
h = hi and h = ^2 be the tangency points (values hi and ^2 are obtained by solving the 
equations <h(/ii) = hi<F'(hi) and <h(h 2 ) = (^2 — l)<h'(h 2 ) -|- 1). According to [I9] solution 
of the Cauchy problem takes the form 

r <F'(hi), he[0,hi) 

e=< <f>'(/i), he[hi,h2] 

[ ‘h'(/i2), h G (^2,1] 


where ^ = {x — xo)/t is the self-similar variable. This means that the solution is presented 
by centred rarefaction wave, which is bounded by two “sonic” shocks. As we show below 
the numerical results of the growth of viscous fingers in the framework of Eqs. dl]), ([5]) are 
in good agreement with the self-similar solution of Eq. flT^ for a suitable choice of the 
parameter k. 
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Time 10 


Time 20 



Figure 11: Comparison of the results for Eqs. (jl]), (E]) (concentration c in two colour is 
presented) and kinematic-wave model flT^ with k = 0.45 (white solid curve). Here we 
choose yUi = 2, /i 2 = 8, and U = 1. 

6.2 Comparing with the Koval model and numerical results 

Let us verify kinematic-wave model flT^ by comparing with the well-approved Koval 
model flThl) . where instead of M effective viscosity ratio Me given by formula flTTjl is used. 
We choose the following parameters: /ii = 2, /i 2 = 8. It means that the viscosity ratio 
M = 4 and according to flT7)) Me = 1.417. Self-similar solution of the Koval model for 
these parameters is shown in Fig. [10] (solid curve). Corresponding solution of the scalar 
conservation law (IT^ (with k = 0.45) is presented in Fig. [10] by dashed curve. As we 
can see in the framework of these models the growth rates of viscous hngers are similar. 
Although the parameter k is a function of M, the values of k{M) are weakly vary for 
M G (1,10). Therefore, for the moderate viscosity ratios M (less than 10) we can assume 
that K = 0.45. Moreover, the proposed model flT^ describes better the structure of the 
huger (its thickness near the tip) than the Koval model flT5|l . 

Further we compare the results of numerical simulations of the formation of viscous 
hngers on the basis of 2D hyperbolic system of equations (|1]), (EJ) with the results obtained 
by using kinematic-wave model flT^ . The calculations are performed in the Hele-Shaw 
cell with sizes L = 20, M = 2 in the coordinate system moving with the average how 
velocity U = 1 with respect to Ox axis. The viscosities of the huids are equal to /ii = 2 
and fi 2 = 8. At the initial time the interface x = Xq = 10 is perturbed as follows 
X = Xq + Ay (exp(—10(|/ — 77/2)^) — 1/2), where Ay = 0.2. For discretization with 
respect to x and y we use 400 and 50 nodes respectively (calculation with hner resolution 
leads to the same result). On the left and right boundaries of the computational domain 
the rehection conditions are imposed. In 2D model (jl]), (|5]) the following dependence for 
viscosity is used: /i(c) = The calculations are carried out for /? = 1, po = 1 and 

the sound velocity Cq = 75. In this case the change in density is not more than 0.15%. At 
the same time the condition = 0 is fulhlled with high accuracy. It means that the above 
proposed ID model is suitable for describing of such flow. Given above perturbation of 
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the interface leads to the formation of a single hnger which is symmetric with respect to 
the line y = \. The resnlts of the concentration c calculations using model (jl]), (jS]) at 
time f = 10 and f = 20 are shown in Fig. [TT] The distribution of c is presented in two 
colours (blue for c < 1/2 and brown for c > 1/2). 

Self-similar solution of kinematic-wave model flT^ with k = 0.45 is shown in Fig. [TT] 
at t = 10 and f = 20 by white solid lines respectively. The curves h + 1 and 1 — h for the 
model are given for a correct comparison with 2D calculations. Here self-similar variable 
^ is replaced hj ^ + U. It corresponds to a transition in a moving coordinate system. The 
hgure shows the propagation velocity and thickness of the viscous hnger obtained by the 
kinematic-wave model agree well with the two-dimensional calculations. 

7 Conclusion 

We derive nonlinear hyperbolic system of equations ([4]), ([5]) describing the how of slightly 
compressible multicomponent huid of diherent viscosity in a Hele-Shaw cell. On the 
basis of these equations simulation of jet how and development of viscous hngers during 
the displacement process are performed. Calculations show that the proposed model 
reproduces the characteristic features of the how associated with the development of 
Kelvin-Helmholtz and Sahman-Taylor instabilities (see Fig. |2l |3] and S]). In the case 
of preferential how in the x-direction the pressure varies slightly in the ^/-direction that 
allows to apply model of long-wave approximation ([6]) and to consider the class of layered 
hows described by Eqs. d?]). Using various simplihcations of the system (linearisation of 
the momentum equations, application of the Darcy’s law) we construct a hierarchy of 
mathematical models of viscosity-stratihed how in a Hele-Shaw cell. These ID models 
are suitable for description of the main features of the two-dimensional how. Stationary 
solutions of Eqs. ([^ obtained for the three-layer how are in good agreement with the 
calculations of the how on the basis of 2D model (jl]), ([1|) (see Fig. Hj). In the framework 
of the three-layer how (systems fll2p and fllbp i the interpretation of the Sahman-Taylor 
instability is given. Solutions of Eqs. (na and flTHl) illustrating the formation of viscous 
hngers are constructed (Fig. [7|and[H|). However, these models fail to predict the correct 
growth rate of the hngers. 

We propose modihcation of the layered how model in order to agree with this be¬ 
haviour. The model is obtained by including friction between the huid layers flTSl) . This 
provides a non-zero thickness of the hngertip and under some assumptions allows one to 
describe the evolution of viscous hngers on the basis of the scalar equation with non-convex 
hux (ITUp . Although the proposed equation flT^ involves empirical parameter, this model 
reveals the physical mechanism to ensure correct propagation velocity and structure of 
viscous hngers. As it can be seen from Fig. [TTjthe velocity of propagation and the thick¬ 
ness of the hngers in proposed model flT^ (if the empirical parameter k is appropriately 
specihed) are in fairly good agreement with calculations based on two-dimensional equa¬ 
tions dlj), dSj)- Comparison with the Koval model (see Fig. ITUl) conhrms the correctness of 
the results. 
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